'''
Created on Aug 17, 2012

@author: tcognata
'''
import math

class calibration(object):
    def __init__(self):
        #Standard set of parameters, perfect calibration
        self.set(alpha=0.003850,beta=0.10863,delta=1.4999,R0=100.0,T0=0.0,Roffset=0.0)
        
    def set(self,**dicty):
        for key,value in dicty.items():
            if (key not in dir(self)) and (key.lower() in [val.lower() for val in dir(self)]):
                key=next(val for val in dir(self) if val.lower()==key.lower())
            self.__setattr__(key,value)
    
    @property
    def alpha(self):
        return self._alpha
    @alpha.setter
    def alpha(self,value):
        self._alpha=value
        
    @property
    def beta(self):
        return self._beta
    @beta.setter
    def beta(self,value):
        self._beta=value
        
    @property
    def delta(self):
        return self._delta
    @delta.setter
    def delta(self,value):
        self._delta=value
        
    @property
    def R0(self):
        return self._R0
    @R0.setter
    def R0(self,value):
        self._R0=value
        
    @property
    def Roffset(self):
        return self._Roffset
    @Roffset.setter
    def Roffset(self,value):
        self._Roffset=value 
        
class RTD(object):
    R_units='ohm'
    T_units='C'
    def __init__(self):
        self.Tout=True
        self.units='C'
        self.cal=calibration()
    
    def __call__(self,R):
        self.R=R
        return self.Temperature(R)
          
    @property
    def A(self):
        return self.cal.alpha*(1+self.cal.delta/100.0)
    @property
    def B(self):
        return -1*self.cal.alpha*self.cal.delta/(100.0*100.0)
    @property
    def C(self):
        return -1*self.cal.alpha*self.cal.beta/(100.0**4.0)
    
    @property
    def R(self):
        if self.Tout:
            return self._resistance
        return self.Resistance(self.T)
    @R.setter
    def R(self,R):
        del self.T
        self._resistance=R
    @R.deleter
    def R(self):
        self.Tout=False

    def Resistance(self,T):
        R0=self.cal.R0
        alpha=self.cal.alpha
        delta=self.cal.delta
        if T>0.0:
            beta=0.0 
        else: 
            beta=self.cal.beta
        S=T/100.0
        return R0+R0*alpha*(T-delta*(S-1.0)*(S)-beta*(S-1.0)*(S)**3.0)
      
    @property
    def T(self):
        if not self.Tout:
            return self._temperature
        return self.Temperature(self.R)
    @T.setter
    def T(self,T):
        del self.R
        self._temperature=T
    @T.deleter
    def T(self):
        self.Tout=True
    
    def Temperature(self,R):
        R0=self.cal.R0
        c=[R0*value for value in [(self.C),(-100*self.C),(self.B),(self.A),(1-R/R0)]]
        #The quadratic solution for resistances greater than R0 (0C)
        if (R >= R0):
            V=((c[3]**2)-(4*c[2]*(R0-R)))
            if (V < 0.0):#raise a warning if sqrt(V) will be complex
                func=lambda T: abs(self.Resistance(T)-R); range=(600,1000)
                warning="Resistance value {} is outside the function range.".format(R)
                return self.Warning(warning,func,range)
            return (-c[3]+math.sqrt(V))/(2*c[2])
        #The quartic solution (per wikipedia) for resistances less than R0 (0C)
        alpha=-(3*c[1]**2)/(8*c[0]**2)+c[2]/c[0]
        beta=(c[1]**3)/(8*c[0]**3)-(c[1]*c[2])/(2*c[0]**2)+(c[3]/c[0])
        gamma=-(3*c[1]**4)/(256*c[0]**4)+(c[2]*c[1]**2)/(16*c[0]**3)-(c[1]*c[3])/(4*c[0]**2)+(c[4]/c[0])
        #if beta==0.0:
        #    simpler solution, see page
        P=-((alpha**2)/12)-gamma
        Q=-((alpha**3)/108)+((alpha*gamma)/3)-((beta**2)/8)
        S=((Q**2)/4)+((P**3)/27)
        if (S < 0.0): #raise a warning if sqrt(S) will be complex
            func=lambda T: abs(self.Resistance(T)-R); range=(-200,-100)
            warning="Resistance value {} is outside the function range.".format(R)
            return self.Warning(warning,func,range)
        R=-(Q/2)+math.sqrt(S) #either sign of the square root will do...
        U=math.pow(R,(1/3))   
        if U==0.0:
            y=-((5/6)*alpha)+U-math.pow(Q,(1/3))
        else:
            y=-((5/6)*alpha)+U-(P/(3*U))   
        W=math.sqrt((alpha+2*y))
        return -(c[1]/(4*c[0]))-W/2+math.sqrt(-(3*alpha+2*y+2*beta/W))/2 
        #For edification:
        #Following is the solution to all parts of the quartic        
        #parts=[-(c[1]/(4*c[0])),W,3*alpha+2*y,2*beta/W]
        #the roots are
        #(import cmath)
        #x1=parts[0]+parts[1]/2+cmath.sqrt(-(parts[2]+parts[3]))/2
        #x2=parts[0]+parts[1]/2-cmath.sqrt(-(parts[2]+parts[3]))/2
        #x3=parts[0]-parts[1]/2+cmath.sqrt(-(parts[2]-parts[3]))/2
        #x4=parts[0]-parts[1]/2-cmath.sqrt(-(parts[2]-parts[3]))/2
        #return x3.real
        
    def Warning(self,message,func=None,range=None):
        import warnings
        warnings.warn(message,)
        if func is not None:
            import scipy.optimize.optimize as op
            return op.brent(func,brack=range)
        
        
def calibrate(R100,R0=100.0,T0=0.0,T100=100.0,Th=None,Rh=None, Tl=None, Rl=None):
    retval={}
    alpha=(R100-R0)/R0/(T100-T0)
    retval['alpha']=alpha
    if (Th and Rh) is not None:
        delta=(Th-((Rh-R0)/(R0*alpha)))/(((Th/100)-1)*(Th/100))
        retval['delta']=delta    
    if (Tl and Rl) is not None:
        beta=(Tl-((Rl-R0)/(R0*alpha)+delta*((Tl/100)-1)*(Tl/100)))/(((Tl/100)-1)*(Tl/100)**3)
        retval['beta']=beta
    
    return retval

if __name__=='__main__':
    import timeit
    item=RTD()

    #item.cal.R0=99.97833
    #item.cal.set(**calibrate(R0=99.97833,T0=-0.01203,R100=138.42492,T100=99.94699,Th=69.96125,Rh=126.97076,Tl=-19.98702,Rl=92.22590))
    #item.cal.set(**calibrate())
    res=[138.42492,126.970076,115.47368,99.97833,96.13693,92.22590]
    a=item(138.42492)
    b=item(92.22590)
    temps=[item(value) for value in res]
    [print(item(value)) for value in res]
    #print('\n\n')
    test1=timeit.Timer('[item(value) for value in res]',"from __main__ import item, res")
    test2=timeit.Timer('[item.Temperature(value) for value in res]',"from __main__ import item, res")
    test3=timeit.Timer('[item.Resistance(value) for value in temps]',"from __main__ import item,temps")
    n=10000
    print('timing for {} calls in the following cases:'.format(n))
    print('object call timing for temperature:    ',test1.timeit(n))
    print('object function timing for temperature:', test2.timeit(n))
    print('object function timing for resistance: ',test3.timeit(n))
    #65 low, 5153 high
    [print(value,item.Temperature(value),item.Resistance(item.Temperature(value))) for value in range(80,350,9)]
    